
\chapter{TAO Solvers}
TAO includes a variety of optimization algorithms for several
classes of problems (unconstrained, bound-constrained,
and PDE-constrained minimization, nonlinear least-squares, and complementarity).
The TAO algorithms for solving these problems are detailed in this section, 
a particular algorithm can chosen by using the \texttt{TaoSetType()} function
or using the command line arguments \texttt{-tao\_type <name>}.
For those interested in extending these algorithms or using new ones, 
please see Chapter~\ref{chapter:addsolver} for more information.


\section{Unconstrained Minimization}
\label{chapter:unconstrained}
Unconstrained minimization is used to minimize a function of many variables
without any constraints on the variables, such as bounds.  The methods 
available in TAO for solving these problems can be classified according
to the amount of derivative information required:
\begin{enumerate}
\item Function evaluation only -- Nelder-Mead method ({\tt tao\_nm})
\item Function and gradient evaluations -- limited-memory, variable-metric 
method ({\tt tao\_lmvm}) and nonlinear conjugate gradient method 
({\tt tao\_cg})
\item Function, gradient, and Hessian evaluations -- Newton line search 
method ({\tt tao\_nls}) and Newton trust-region method ({\tt tao\_ntr})
\end{enumerate}
The best method to use depends on the particular problem being solved
and the accuracy required in the solution.  If a Hessian evaluation 
routine is available, then the Newton line search and Newton trust-region 
methods will likely perform best.  When a Hessian evaluation routine
is not available, then the limited-memory, variable-metric method is 
likely to perform best.  The Nelder-Mead method should be used only
as a last resort when no gradient information is available.

Each solver has a set of options associated with it that can be set with 
command line arguments.  These algorithms and the 
associated options are briefly discussed in this chapter.

\subsection{Nelder-Mead Method}
The Nelder-Mead algorithm \cite{nelder.mead:simplex} is a direct search method for finding a local
minimum of a function $f(x)$.  This algorithm does not require any gradient or Hessian 
information of $f$ and therefore has some expected advantages and disadvantages compared
to the other TAO solvers.  The obvious advantage is that it is easier to write an 
application when no derivatives need to be calculated.  The downside is that this algorithm can
be slow to converge or can even stagnate, and it performs poorly for large numbers of variables.

This solver keeps a set of $N+1$ sorted vectors ${x_1,x_2,\ldots,x_{N+1}}$ and their corresponding 
objective function values $f_1 \leq f_2 \leq \ldots \leq f_{N+1}$.  At each iteration, $x_{N+1}$ is removed from
the set and replaced with 
\[
x(\mu) = (1+\mu) \frac{1}{N} \sum_{i=1}^N x_i - \mu x_{N+1},
\]
where $\mu$ can be one of
${\mu_0,2\mu_0,\frac{1}{2}\mu_0,-\frac{1}{2}\mu_0}$ depending on the
values of 
each possible $f(x(\mu))$.

The algorithm terminates when the residual  $f_{N+1} - f_1$ becomes sufficiently small.  Because of 
the way new vectors can be added to the sorted set, 
the minimum function value and/or the residual may not be impacted at each iteration.

Two options can be set specifically for the Nelder-Mead algorithm:
\begin{description}
\item[\tt -tao\_nm\_lamda <value>] sets the initial set of vectors ($x_0$ plus 
{\tt value} in each coordinate direction); the default value is $1$.  
\item[\tt -tao\_nm\_mu <value>] sets the value of $\mu_0$; the default 
is $\mu_0=1$.
\end{description}

\subsection{Limited-Memory, Variable-Metric Method}\sindex{line search}\sindex{gradients}

The limited-memory, variable-metric method computes a positive definite
approximation to the Hessian matrix from a limited number of previous
iterates and gradient evaluations.  A direction is then obtained by
solving the system of equations
\[
H_k d_k = -\nabla f(x_k),
\]
where $H_k$ is the Hessian approximation obtained by using the BFGS update 
formula.  The inverse of $H_k$ can readily be applied to obtain the 
direction $d_k$.  Having obtained the direction, a Mor\'{e}-Thuente 
line search is applied to compute a step length, $\tau_k$, that 
approximately solves the one-dimensional optimization problem
\[
\min_\tau f(x_k + \tau d_k).
\]
The current iterate and Hessian approximation are updated, and the process
is repeated until the method converges.  This algorithm is the default 
unconstrained minimization solver and can be selected by using the 
TAO solver {\tt tao\_lmvm}.  For best efficiency, function and gradient 
evaluations should be performed simultaneously when using this algorithm.

The primary factors determining the behavior of this algorithm are the type 
of Hessian approximation used, the number of vectors stored for the approximation 
and the initialization/scaling of the approximation. These options can be 
configured using the {\tt -tao\_lmvm\_mat\_lmvm} prefix. For further detail, 
we refer the reader to the {\tt MATLMVM} matrix type definitions in the PETSc 
Manual.

The LMVM algorithm also allows the user to define a custom initial 
Hessian matrix $H_{0,k}$ through the interface function {\tt TaoLMVMSetH0()}. 
This user-provided initialization overrides any other scalar or diagonal 
initialization inherent to the LMVM approximation. The provided $H_{0,k}$ must 
be a PETSc {\tt Mat} type object that represents a positive-definite matrix. 
The approximation prefers {\tt MatSolve()} if the provided matrix has 
{\tt MATOP\_SOLVE} implemented. Otherwise, {\tt MatMult()} is used in a KSP 
solve to perform the inversion of the user-provided initial Hessian.

In applications where {\tt TaoSolve()} on the LMVM algorithm is repeatedly called 
to solve similar or related problems, {\tt -tao\_lmvm\_recycle} flag can be used 
to prevent resetting the LMVM approximation between subsequent solutions. This 
recycling also avoids one extra function and gradient evalution, instead re-using 
the values already computed at the end of the previous solution.

This algorithm will be deprecated in the next version and replaced by the 
bounded quasi-Newton Line Search (BQNLS) algorithm that can solve both bound constrained 
and unconstrained problems.

\subsection{Nonlinear Conjugate Gradient Method}\sindex{line search}\sindex{gradients}

The nonlinear conjugate gradient method can be viewed as an extension of the 
conjugate gradient method for solving symmetric, positive-definite linear 
systems of equations.  This algorithm requires only function and gradient 
evaluations as well as a line search.  The TAO implementation uses a 
Mor\'{e}-Thuente line search to obtain the step length.  The nonlinear 
conjugate gradient method can be selected by using the TAO solver 
{\tt tao\_cg}.  For the best efficiency, function and gradient evaluations 
should be performed simultaneously when using this algorithm.

Five variations are currently supported by the TAO implementation: the 
Fletcher-Reeves method, the Polak-Ribi\'ere method, the Polak-Ribi\'ere-Plus 
method~\cite{NW99}, the Hestenes-Stiefel method, and the Dai-Yuan method.  
These conjugate gradient methods can be specified by using the command line 
argument {\tt -tao\_cg\_type <fr,pr,prp,hs,dy>}, respectively.  The default 
value is {\tt prp}.  

The conjugate gradient method incorporates automatic restarts when successive 
gradients are not sufficiently orthogonal.  TAO measures the orthogonality by 
dividing the inner product of the gradient at the current point and the 
gradient at the previous point by the square of the Euclidean norm of 
the gradient at the current point.  When the absolute value of this 
ratio is greater than $\eta$, the algorithm restarts using the gradient 
direction.  The parameter $\eta$ can be set by using the command line argument 
{\tt -tao\_cg\_eta <real>}; 0.1 is the default value.  

\subsection{Newton Line Search Method}\sindex{Newtons method}\sindex{line search}

The Newton line search method solves the symmetric system of equations
\[
H_k d_k = -g_k
\]
to obtain a step $d_k$, where $H_k$ is the Hessian of the objective function
at $x_k$ and $g_k$ is the gradient of the objective function at $x_k$.
For problems where the Hessian matrix is indefinite, the perturbed system
of equations
\[
(H_k + \rho_k I) d_k = -g_k
\]
is solved to obtain the direction, where $\rho_k$ is a positive constant.
If the direction computed is not a descent direction, the (scaled) steepest 
descent direction is used instead.  Having obtained the direction, 
a Mor\'{e}-Thuente line search is applied to obtain a step length, 
$\tau_k$, that approximately solves the one-dimensional optimization 
problem
\[
\min_\tau f(x_k + \tau d_k).
\]
The Newton line search method can be selected by using the TAO 
solver {\tt tao\_nls}.
The options available for this solver are listed in
Table~\ref{table:nlsoptions}.
For the best efficiency, function and gradient evaluations should be 
performed simultaneously when using this algorithm.

\begin{table}[h]
\caption{Summary of {\tt nls} options}
\label{table:nlsoptions}
\vspace{5mm}
\begin{tabular}{l|p{1.5in}|l|p{2.0in}}
\hline
Name \rule{0pt}{2.6ex} & Value & Default & Description \\
\hline
{\tt -tao\_nls\_ksp\_type} \rule{0pt}{2.6ex} & cg, nash, stcg, gltr, gmres, ... & stcg & Type of Krylov subspace method to use when solving linear system \\
{\tt -tao\_nls\_pc\_type} & none, jacobi, icc, ilu, lmvm & lmvm & Type of preconditioner to use when solving linear system \\
{\tt -tao\_nls\_sval} & real & $0$ & Initial perturbation value \\
{\tt -tao\_nls\_imin} & real & $10^{-4}$ & Minimum initial perturbation value \\
{\tt -tao\_nls\_imax} & real & $100$ & Maximum initial perturbation value \\
{\tt -tao\_nls\_imfac} & real & $0.1$ & Factor applied to norm of gradient when initializing perturbation \\
{\tt -tao\_nls\_pmax} & real & $100$ & Maximum perturbation when increasing value \\
{\tt -tao\_nls\_pgfac} & real & $10$ & Growth factor applied to perturbation when increasing value \\
{\tt -tao\_nls\_pmgfac} & real & $0.1$ & Factor applied to norm of gradient when increasing perturbation \\
{\tt -tao\_nls\_pmin} & real & $10^{-12}$ & Minimum perturbation when decreasing value; smaller values set to zero \\
{\tt -tao\_nls\_psfac} & real & $0.4$ & Shrink factor applied to perturbation when decreasing value \\
{\tt -tao\_nls\_pmsfac} & real & $0.1$ & Factor applied to norm of gradient when decreasing perturbation \\
{\tt -tao\_nls\_init\_type} & constant, direction, interpolation &
interpolation & Method used to initialize trust-region radius when using
{\tt nash}, {\tt stcg}, or {\tt gltr} \\
\hline
\end{tabular}
\end{table}
\afterpage{\clearpage}

\begin{table}[h]
\caption{Summary of {\tt nls} options (continued)}
\vspace{5mm}
\begin{tabular}{l|p{1.5in}|l|p{2.0in}}
\hline
Name \rule{0pt}{2.6ex} & Value & Default & Description \\
\hline
{\tt -tao\_nls\_mu1\_i} \rule{0pt}{2.6ex} & real & 0.35 & $\mu_1$ in {\tt interpolation} init \\
{\tt -tao\_nls\_mu2\_i} & real & 0.50 & $\mu_2$ in {\tt interpolation} init \\
{\tt -tao\_nls\_gamma1\_i} & real & 0.0625 & $\gamma_1$ in {\tt interpolation} init \\
{\tt -tao\_nls\_gamma2\_i} & real & 0.50 & $\gamma_2$ in {\tt interpolation} init \\
{\tt -tao\_nls\_gamma3\_i} & real & 2.00 & $\gamma_3$ in {\tt interpolation} init \\
{\tt -tao\_nls\_gamma4\_i} & real & 5.00 & $\gamma_4$ in {\tt interpolation} init \\
{\tt -tao\_nls\_theta\_i} & real & 0.25 & $\theta$ in {\tt interpolation} init \\
{\tt -tao\_nls\_update\_type} & step, reduction, interpolation & step &
Method used to update trust-region radius when using {\tt nash}, {\tt
stcg}, or {\tt gltr} \\
{\tt -tao\_nls\_nu1} & real & 0.25 & $\nu_1$ in {\tt step} update \\
{\tt -tao\_nls\_nu2} & real & 0.50 & $\nu_2$ in {\tt step} update \\
{\tt -tao\_nls\_nu3} & real & 1.00 & $\nu_3$ in {\tt step} update \\
{\tt -tao\_nls\_nu4} & real & 1.25 & $\nu_4$ in {\tt step} update \\
{\tt -tao\_nls\_omega1} & real & 0.25 & $\omega_1$ in {\tt step} update \\
{\tt -tao\_nls\_omega2} & real & 0.50 & $\omega_2$ in {\tt step} update \\
{\tt -tao\_nls\_omega3} & real & 1.00 & $\omega_3$ in {\tt step} update \\
{\tt -tao\_nls\_omega4} & real & 2.00 & $\omega_4$ in {\tt step} update \\
{\tt -tao\_nls\_omega5} & real & 4.00 & $\omega_5$ in {\tt step} update \\
{\tt -tao\_nls\_eta1} & real & $10^{-4}$ & $\eta_1$ in {\tt reduction} update \\
{\tt -tao\_nls\_eta2} & real & 0.25 & $\eta_2$ in {\tt reduction} update \\
{\tt -tao\_nls\_eta3} & real & 0.50 & $\eta_3$ in {\tt reduction} update \\
{\tt -tao\_nls\_eta4} & real & 0.90 & $\eta_4$ in {\tt reduction} update \\
{\tt -tao\_nls\_alpha1} & real & 0.25 & $\alpha_1$ in {\tt reduction} update \\
{\tt -tao\_nls\_alpha2} & real & 0.50 & $\alpha_2$ in {\tt reduction} update \\
{\tt -tao\_nls\_alpha3} & real & 1.00 & $\alpha_3$ in {\tt reduction} update \\
{\tt -tao\_nls\_alpha4} & real & 2.00 & $\alpha_4$ in {\tt reduction} update \\
{\tt -tao\_nls\_alpha5} & real & 4.00 & $\alpha_5$ in {\tt reduction} update \\
{\tt -tao\_nls\_mu1} & real & 0.10 & $\mu_1$ in {\tt interpolation} update \\
{\tt -tao\_nls\_mu2} & real & 0.50 & $\mu_2$ in {\tt interpolation} update \\
{\tt -tao\_nls\_gamma1} & real & 0.25 & $\gamma_1$ in {\tt interpolation} update \\
{\tt -tao\_nls\_gamma2} & real & 0.50 & $\gamma_2$ in {\tt interpolation} update \\
{\tt -tao\_nls\_gamma3} & real & 2.00 & $\gamma_3$ in {\tt interpolation} update \\
{\tt -tao\_nls\_gamma4} & real & 4.00 & $\gamma_4$ in {\tt interpolation} update \\
{\tt -tao\_nls\_theta} & real & 0.05 & $\theta$ in {\tt interpolation} update \\
\hline
\end{tabular}
\end{table}
\afterpage{\clearpage}

The system of equations is approximately solved by applying the conjugate 
gradient method, Nash conjugate gradient method, Steihaug-Toint conjugate 
gradient method, generalized 
Lanczos method, or an alternative Krylov subspace method 
supplied by PETSc.  The method used to solve the systems of equations is 
specified with the command line argument 
{\tt -tao\_nls\_ksp\_type <cg,nash,stcg,gltr,gmres,...>}; {\tt stcg} 
is the default. See the PETSc manual for further information on changing 
the behavior of the linear system solvers.  

A good preconditioner reduces the number of iterations required to
solve the linear system of equations.  For the conjugate gradient
methods and generalized Lanczos method, this preconditioner must be
symmetric and positive definite.  The available options are to use no
preconditioner, the absolute value of the diagonal of the Hessian
matrix, a limited-memory BFGS approximation to the Hessian matrix, or
one of the other preconditioners provided by the PETSc package.  These
preconditioners are specified by the command line arguments 
{\tt -tao\_nls\_pc\_type <none,jacobi,icc,ilu,lmvm>}, respectively. The
default is the {\tt lmvm} preconditioner, which uses a BFGS approximation 
of the inverse Hessian. See the PETSc manual for further information on 
changing the behavior of the preconditioners.

The perturbation $\rho_k$ is added when the direction returned by the
Krylov subspace method is not a descent direction, the Krylov method
diverged due to an indefinite preconditioner or matrix, or a direction of 
negative curvature was found.  In the last two cases, if the step returned
is a descent direction, it is used during the line search.  Otherwise, a
steepest descent direction is used during the line search.  The perturbation
is decreased as long as the Krylov subspace method reports success and 
increased if further problems are encountered.  There are three cases:
initializing, increasing, and decreasing the perturbation.  These cases
are described below.
\begin{enumerate}
\item If $\rho_k$ is zero and a problem was detected with either the
direction or the Krylov subspace method, the perturbation is initialized to
\[
\rho_{k+1} = \mbox{median}\left\{\mbox{imin}, \mbox{imfac} * \|g(x_k)\|, \mbox{imax}\right\},
\]
where $g(x_k)$ is the gradient of the objective function and {\tt imin} is set 
with the command line argument {\tt -tao\_nls\_imin <real>} with a default 
value of $10^{-4}$, {\tt imfac} by {\tt -tao\_nls\_imfac} with a default 
value of 0.1, and {\tt imax} by {\tt -tao\_nls\_imax} with a default value 
of 100.  
When using the {\tt gltr} method to solve the system of equations, an
estimate of the minimum eigenvalue $\lambda_1$ of the Hessian matrix 
is available.  This value is used to initialize the perturbation to
$\rho_{k+1} = \max\left\{\rho_{k+1}, -\lambda_1\right\}$ in this
case.
\item If $\rho_k$ is nonzero and a problem was detected with either the 
direction or Krylov subspace method, the perturbation is increased to 
\[
\rho_{k+1} = \min\left\{\mbox{pmax}, \max\left\{\mbox{pgfac} * \rho_k, \mbox{pmgfac} * \|g(x_k)\|\right\}\right\},
\]
where $g(x_k)$ is the gradient of the objective function and {\tt pgfac} is set 
with the command line argument {\tt -tao\_nls\_pgfac}
with a default value of 10, {\tt pmgfac} by {\tt -tao\_nls\_pmgfac} with a
default value of 0.1, and {\tt pmax} by {\tt -tao\_nls\_pmax} with a default
value of 100.
\item If $\rho_k$ is nonzero and no problems were detected with either
the direction or Krylov subspace method, the perturbation is decreased to
\[
\rho_{k+1} = \min\left\{\mbox{psfac} * \rho_k, \mbox{pmsfac} * \|g(x_k)\|\right\},
\]
where $g(x_k)$ is the gradient of the objective function, {\tt psfac} is set 
with the command line argument {\tt -tao\_nls\_psfac}
with a default value of 0.4, and {\tt pmsfac} is set by {\tt -tao\_nls\_pmsfac}
with a default value of 0.1.  Moreover, if $\rho_{k+1} < \mbox{pmin}$, then 
$\rho_{k+1} = 0$, where {\tt pmin} is set with the command line argument 
{\tt -tao\_nls\_pmin} and has a default value of $10^{-12}$.
\end{enumerate}
Near a local minimizer to the unconstrained optimization problem, the
Hessian matrix will be positive-semidefinite; the perturbation will 
shrink toward zero, and one would eventually observe a superlinear
convergence rate.

When using {\tt nash}, {\tt stcg}, or {\tt gltr} to solve the linear systems 
of equation,
a trust-region radius needs to be initialized and updated.  This
trust-region
radius simultaneously limits the size of the step computed and reduces the 
number of iterations of the conjugate gradient method.  The method for 
initializing the trust-region radius is set with the command line argument 
{\tt -tao\_nls\_init\_type <constant,direction,interpolation>};
{\tt interpolation}, which chooses an initial value based on the 
interpolation scheme found in \cite{CGT}, is the default.  This
scheme performs a number of function and gradient evaluations to determine 
a radius such that the reduction predicted by the quadratic model along the 
gradient direction coincides with the actual reduction in the nonlinear 
function.  The iterate obtaining the best objective function value is 
used as the starting point for the main line search algorithm.  The 
{\tt constant} method initializes the trust-region radius by using 
the value specified with the {\tt -tao\_trust0 <real>} command line 
argument, where the default value is 100.  The {\tt direction} technique 
solves the first quadratic optimization problem by using a standard 
conjugate gradient method and initializes the trust region to 
$\|s_0\|$.

The method for updating the trust-region radius is set with the 
command line argument 
{\tt -tao\_nls\_update\_type <step,reduction,interpolation>}; {\tt step} 
is the default.  The {\tt step} method updates the trust-region 
radius based on the value of $\tau_k$.  In particular,
\[
\Delta_{k+1} = \left\{\begin{array}{ll}
\omega_1 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \tau_k \in [0, \nu_1) \\
\omega_2 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \tau_k \in [\nu_1, \nu_2) \\
\omega_3 \Delta_k & \mbox{if } \tau_k \in [\nu_2, \nu_3) \\
\mbox{max}(\Delta_k, \omega_4 \|d_k\|) & \mbox{if } \tau_k \in [\nu_3, \nu_4) \\
\mbox{max}(\Delta_k, \omega_5 \|d_k\|) & \mbox{if } \tau_k \in [\nu_4, \infty),
\end{array}
\right.
\]
where $0 < \omega_1 < \omega_2 < \omega_3 = 1 < \omega_4 < \omega_5$ and
$0 < \nu_1 < \nu_2 < \nu_3 < \nu_4$ are constants.  The {\tt reduction} 
method computes the ratio of the actual reduction in the objective function 
to the reduction predicted by the quadratic model for the full step, 
$\kappa_k = \frac{f(x_k) - f(x_k + d_k)}{q(x_k) - q(x_k + d_k)}$, where 
$q_k$ is the quadratic model.  The radius is then updated as
\[
\Delta_{k+1} = \left\{\begin{array}{ll}
\alpha_1 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \kappa_k \in (-\infty, \eta_1) \\
\alpha_2 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \kappa_k \in [\eta_1, \eta_2) \\
\alpha_3 \Delta_k & \mbox{if } \kappa_k \in [\eta_2, \eta_3) \\
\mbox{max}(\Delta_k, \alpha_4 \|d_k\|) & \mbox{if } \kappa_k \in [\eta_3, \eta_4) \\
\mbox{max}(\Delta_k, \alpha_5 \|d_k\|) & \mbox{if } \kappa_k \in [\eta_4, \infty),
\end{array}
\right.
\]
where $0 < \alpha_1 < \alpha_2 < \alpha_3 = 1 < \alpha_4 < \alpha_5$ and
$0 < \eta_1 < \eta_2 < \eta_3 < \eta_4$ are constants.  The {\tt interpolation}
method uses the same interpolation mechanism as in the initialization to
compute a new value for the trust-region radius.

This algorithm will be deprecated in the next version and replaced by the 
bounded Newton Line Search (BNLS) algorithm that can solve both bound constrained 
and unconstrained problems.

\subsection{Newton Trust-Region Method}\sindex{Newton method}\sindex{trust
region}

The Newton trust-region method solves the constrained quadratic programming
problem
\[
\begin{array}{ll}
\min_d  & \frac{1}{2}d^T H_k d  + g_k^T d \\
\mbox{subject to} & \|d\| \leq \Delta_k
\end{array}
\]
to obtain a direction $d_k$, where $H_k$ is the Hessian of the objective 
function at $x_k$, $g_k$ is the gradient of the objective function at $x_k$, 
and $\Delta_k$ is the trust-region radius.  If $x_k + d_k$ sufficiently 
reduces the nonlinear objective function, then the step is accepted, and the 
trust-region radius is updated.  However, if $x_k + d_k$ does not
sufficiently
reduce the nonlinear objective function, then the step is rejected, the 
trust-region radius is reduced, and the quadratic program is re-solved 
by using the updated trust-region radius. The Newton trust-region method 
can be set by using the TAO solver {\tt tao\_ntr}. The options available 
for this solver are listed in Table~\ref{table:ntroptions}. For the best 
efficiency, 
function and gradient evaluations should be performed separately when 
using this algorithm.

\begin{table}[h]
\caption{Summary of {\tt ntr} options}
\label{table:ntroptions}
\vspace{5mm}
\begin{tabular}{l|p{1.5in}|l|p{2.0in}}
\hline
Name \rule{0pt}{2.6ex} & Value & Default & Description \\
\hline
{\tt -tao\_ntr\_ksp\_type} \rule{0pt}{2.6ex} & nash, stcg, gltr & stcg & Type of Krylov subspace method to use when solving linear system \\
{\tt -tao\_ntr\_pc\_type} & none, jacobi, icc, ilu, lmvm & lmvm & Type of preconditioner to use when solving linear system \\
{\tt -tao\_ntr\_init\_type} & constant, direction, interpolation &
interpolation & Method used to initialize trust-region radius \\
{\tt -tao\_ntr\_mu1\_i} & real & 0.35 & $\mu_1$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_mu2\_i} & real & 0.50 & $\mu_2$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_gamma1\_i} & real & 0.0625 & $\gamma_1$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_gamma2\_i} & real & 0.50 & $\gamma_2$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_gamma3\_i} & real & 2.00 & $\gamma_3$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_gamma4\_i} & real & 5.00 & $\gamma_4$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_theta\_i} & real & 0.25 & $\theta$ in {\tt interpolation} init \\
{\tt -tao\_ntr\_update\_type} & reduction, interpolation & reduction &
Method used to update trust-region radius \\
{\tt -tao\_ntr\_eta1} & real & $10^{-4}$ & $\eta_1$ in {\tt reduction} update \\
{\tt -tao\_ntr\_eta2} & real & 0.25 & $\eta_2$ in {\tt reduction} update \\
{\tt -tao\_ntr\_eta3} & real & 0.50 & $\eta_3$ in {\tt reduction} update \\
{\tt -tao\_ntr\_eta4} & real & 0.90 & $\eta_4$ in {\tt reduction} update \\
{\tt -tao\_ntr\_alpha1} & real & 0.25 & $\alpha_1$ in {\tt reduction} update \\
{\tt -tao\_ntr\_alpha2} & real & 0.50 & $\alpha_2$ in {\tt reduction} update \\
{\tt -tao\_ntr\_alpha3} & real & 1.00 & $\alpha_3$ in {\tt reduction} update \\
{\tt -tao\_ntr\_alpha4} & real & 2.00 & $\alpha_4$ in {\tt reduction} update \\
{\tt -tao\_ntr\_alpha5} & real & 4.00 & $\alpha_5$ in {\tt reduction} update \\
{\tt -tao\_ntr\_mu1} & real & 0.10 & $\mu_1$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_mu2} & real & 0.50 & $\mu_2$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_gamma1} & real & 0.25 & $\gamma_1$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_gamma2} & real & 0.50 & $\gamma_2$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_gamma3} & real & 2.00 & $\gamma_3$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_gamma4} & real & 4.00 & $\gamma_4$ in {\tt interpolation} update \\
{\tt -tao\_ntr\_theta} & real & 0.05 & $\theta$ in {\tt interpolation} update \\
\hline
\end{tabular}
\end{table}
\afterpage{\clearpage}

The quadratic optimization problem is approximately solved by applying
the Nash or Steihaug-Toint conjugate gradient methods or the 
generalized Lanczos 
method to the symmetric system of equations $H_k d = -g_k$.  The method 
used to solve the system of equations is specified with the command line
argument {\tt -tao\_ntr\_ksp\_type <nash,stcg,gltr>}; {\tt stcg} is the default.  
See the PETSc manual for further information on changing the behavior of 
these linear system solvers.  

A good preconditioner reduces the number of iterations required to
compute the direction.  For the Nash and Steihaug-Toint conjugate gradient
methods and generalized Lanczos method, this preconditioner must be
symmetric and positive definite.  The available options are to use no
preconditioner, the absolute value of the diagonal of the Hessian
matrix, a limited-memory BFGS approximation to the Hessian matrix, or
one of the other preconditioners provided by the PETSc package.  These
preconditioners are specified by the command line argument 
{\tt -tao\_ntr\_pc\_type <none,jacobi,icc,ilu,lmvm>}, respectively.  The
default is the {\tt lmvm} preconditioner. See the PETSc manual 
for further information on changing the behavior of the preconditioners.

The method for computing an initial trust-region radius is set with the 
command line arguments {\tt -tao\_ntr\_init\_type <constant,direction,interpolation>};
{\tt interpolation}, which chooses an initial value based on the 
interpolation scheme found in \cite{CGT}, is the default.  This
scheme performs a number of function and gradient evaluations to determine 
a radius such that the reduction predicted by the quadratic model along the 
gradient direction coincides with the actual reduction in the nonlinear 
function.  The iterate obtaining the best objective function value is 
used as the starting point for the main trust-region algorithm.  The 
{\tt constant} method initializes the trust-region radius by using 
the value specified with the {\tt -tao\_trust0 <real>} command line 
argument, where the default value is 100.  The {\tt direction} technique 
solves the first quadratic optimization problem by using a standard 
conjugate gradient method and initializes the trust region to 
$\|s_0\|$.

The method for updating the trust-region radius is set with the 
command line arguments 
{\tt -tao\_ntr\_update\_type <reduction,interpolation>}; {\tt reduction} 
is the default.  The {\tt reduction} method computes the ratio of the 
actual reduction in the objective function to the reduction predicted 
by the quadratic model for the full step, 
$\kappa_k = \frac{f(x_k) - f(x_k + d_k)}{q(x_k) - q(x_k + d_k)}$, where 
$q_k$ is the quadratic model.  The radius is then updated as
\[
\Delta_{k+1} = \left\{\begin{array}{ll}
\alpha_1 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \kappa_k \in (-\infty, \eta_1) \\
\alpha_2 \mbox{min}(\Delta_k, \|d_k\|) & \mbox{if } \kappa_k \in [\eta_1, \eta_2) \\
\alpha_3 \Delta_k & \mbox{if } \kappa_k \in [\eta_2, \eta_3) \\
\mbox{max}(\Delta_k, \alpha_4 \|d_k\|) & \mbox{if } \kappa_k \in [\eta_3, \eta_4) \\
\mbox{max}(\Delta_k, \alpha_5 \|d_k\|) & \mbox{if } \kappa_k \in [\eta_4, \infty),
\end{array}
\right.
\]
where $0 < \alpha_1 < \alpha_2 < \alpha_3 = 1 < \alpha_4 < \alpha_5$ and
$0 < \eta_1 < \eta_2 < \eta_3 < \eta_4$ are constants.  The {\tt interpolation}
method uses the same interpolation mechanism as in the initialization to
compute a new value for the trust-region radius.

This algorithm will be deprecated in the next version and replaced by the 
bounded Newton Trust Region (BNTR) algorithm that can solve both bound constrained 
and unconstrained problems.

\subsection{BMRM}\sindex{BMRM}
The Bundle Method for Regularized Risk Minimization (BMRM)\cite{brmrm} is a 
numerical approach to optimizing an unconstrained objective in the form of 
$f(x) + 0.5 * \lambda \| x \|^2$.  Here $f$ is a convex function that is finite on 
the whole space.  $\lambda$ is a positive 
weight parameter, and $\| x \|$ is the Euclidean norm of $x$.  The algorithm only 
requires a routine which, given an $x$, returns the value of $f(x)$ and the 
gradient of $f$ at $x$.

\subsection{OWL-QN}\sindex{OWLQN}
The Orthant-Wise Limited-memory Quasi-Newton algorithm (OWL-QN)\cite{owlqn}
is a numerical approach to optimizing an unconstrained objective in the 
form of $f(x) + \lambda \|x\|_1$.  Here f is a convex and differentiable 
function, $\lambda$ is a positive weight parameter, and $\| x \|_1$ is the 
$L1$ norm of $x$: $\sum_i |x_i|$.  The algorithm only requires evaluating the 
value of $f$ and its gradient.

\section{Bound-Constrained Optimization}\sindex{bounds}
\label{chapter:bound}

Bound-constrained optimization algorithms solve optimization problems of the form
\[
\begin{array}{ll} \displaystyle
\min_{x} & f(x) \\
\mbox{subject to} & l \leq x \leq u.
\end{array}
\]
These solvers use the bounds on the variables as well as objective
function, gradient, and possibly Hessian information.

For any unbounded variables, the bound value for the associated index can be set 
to {\tt PETSC\_INFINITY} for the upper bound and {\tt PETSC\_NINFINITY} for the lower 
bound. If all bounds are set to infinity, then the bounded algorithms are equivalent 
to their unconstrained counterparts.

Before introducing specific methods, we will first define two projection operations 
used by all bound constrained algorithms.
\begin{itemize}
\item Gradient projection:
\[
\mathfrak{P}(g) = \left\{\begin{array}{ll}
0 & \mbox{if} \; (x \leq l_i \land g_i > 0) \lor (x \geq u_i \land g_i < 0) \\
g_i & \mbox{otherwise}
\end{array}
\right.
\]
\item Bound projection:
\[
\mathfrak{B}(x) = \left\{\begin{array}{ll}
l_i & \mbox{if} \; x_i < l_i \\
u_i & \mbox{if} \; x_i > u_i \\
x_i & \mbox{otherwise}
\end{array}
\right.
\]
\end{itemize}

\subsection{Bounded Newton-Krylov Methods}\sindex{projected Newton's
method}\sindex{trust region}\label{sec:bnk}

TAO features three bounded Newton-Krylov (BNK) class of algorithms, separated by their 
globalization methods: projected line search (BNLS), trust region (BNTR), and trust region 
with a projected line search fall-back (BNTL). They are available via the TAO 
solvers {\tt TAOBNLS}, {\tt TAOBNTR} and {\tt TAOBNTL}, respectively, or the 
{\tt -tao\_type} {\tt bnls}/{\tt bntr}/{\tt bntl} flag.

The BNK class of methods use an active-set approach to solve the symmetric system of equations,
\[
H_k p_k = -g_k,
\]
only for inactive variables in the interior of the bounds. The active-set estimation is based on 
Bertsekas \cite{bertsekas:projected} with the following variable index categories:
\[
\begin{array}{rlll} \displaystyle
\mbox{lower bounded}: & \mathcal{L}(x) & = & \{ i \; : \; x_i \leq l_i + \epsilon \; \land \; g(x)_i > 0 \}, \\
\mbox{upper bounded}: & \mathcal{U}(x) & = & \{ i \; : \; x_i \geq u_i + \epsilon \; \land \; g(x)_i < 0 \}, \\
\mbox{fixed}: & \mathcal{F}(x) & = & \{ i \; : \; l_i = u_i \}, \\
\mbox{active-set}: & \mathcal{A}(x) & = & \{ \mathcal{L}(x) \; \bigcup \; \mathcal{U}(x) \; \bigcup \; \mathcal{F}(x) \}, \\
\mbox{inactive-set}: & \mathcal{I}(x) & = & \{ 1,2,\ldots,n \} \; \backslash \; \mathcal{A}(x).
\end{array}
\]
At each iteration, the bound tolerance is estimated as $\epsilon_{k+1} = \mbox{min}(\epsilon_k, ||w_k||_2)$ 
with $w_k = x_k - \mathfrak{B}(x_k - \beta D_k g_k)$, where the diagonal matrix $D_k$ is an approximation 
of the Hessian inverse $H_k^{-1}$. The initial bound tolerance $\epsilon_0$ and the step length $\beta$ 
have default values of $0.001$ and can be adjusted using {\tt -tao\_bnk\_as\_tol} and 
{\tt -tao\_bnk\_as\_step} flags, respectively. The active-set estimation can be disabled using the option 
{\tt -tao\_bnk\_as\_type none}, in which case the algorithm simply uses the current iterate with no bound 
tolerances to determine which variables are actively bounded and which are free.

BNK algorithms invert the reduced Hessian using a Krylov iterative method. Trust-region conjugate gradient 
methods ({\tt KSPNASH}, {\tt KSPSTCG}, and {\tt KSPGLTR}) are required for the BNTR and BNTL algorithms, 
and recommended for the BNLS algorithm. The preconditioner type can be changed using the 
{\tt -tao\_bnk\_pc\_type} {\tt none}/{\tt ilu}/{\tt icc}/{\tt jacobi}/{\tt lmvm}. The {\tt lmvm} option, which is 
also the default, preconditions the Krylov solution with a {\tt MATLMVM} matrix. The remaining supported 
preconditioner types are default PETSc types. If Jacobi is selected, the diagonal values are safeguarded to be positive. 
{\tt icc} and {\tt ilu} options produce good results for problems with dense Hessians. The LMVM and Jacobi 
preconditioners are also used as the approximate inverse-Hessian in the active-set estimation. 
If neither are available, or if the Hessian matrix does not have {\tt MATOP\_GET\_DIAGONAL} defined, then the 
active-set estimation falls back onto using an identity matrix in place of $D_k$ (this is equivalent to 
estimating the active-set using a gradient descent step).

A special option is available to ``accelerate'' the convergence of the BNK algorithms by taking a finite 
number of BNCG iterations at each Newton iteration. By default, the number of BNCG iterations is set to zero 
and the algorithms do not take any BNCG steps. This can be changed using the option flag 
{\tt -tao\_bnk\_max\_cg\_its <i>}. While this reduces the number of Newton iterations, in practice it simply 
trades off the Hessian evaluations in the BNK solver for more function and gradient evaluations in the BNCG 
solver. However, it may be useful for certain types of problems where the Hessian evaluation is 
disproportionately more expensive than the objective function or its gradient.

\subsubsection{Bounded Newton Line Search (BNLS)}\label{sec:bnls}

BNLS safeguards the Newton step by falling back onto a BFGS, scaled gradient, or gradient 
steps based on descent direction verifications. For problems with indefinite Hessian matrices, 
the step direction is calculated using a perturbed system of equations,
\[
(H_k + \rho_k I)p_k = -g_k,
\]
where $\rho_k$ is a dynamically adjusted positive constant. The step is globalized using a 
projected More-Thuente line search. If a trust-region conjugate gradient method is used for the 
Hessian inversion, the trust radius is modified based on the line search step length.

\subsubsection{Bounded Newton Trust Region (BNTR)}\label{sec:bntr}

BNTR globalizes the Newton step using a trust region method based on the predicted versus 
actual reduction in the cost function. The trust radius is increased only if the accepted 
step is at the trust region boundary. The reduction check features a safeguard for numerical 
values below machine epsilon, scaled by the latest function value, where the full Newton step 
is accepted without modification.

\subsubsection{Bounded Newton Trust Region with Line Search Fall-back (BNTL)}\label{sec:bntr}

BNTL safeguards the trust-region globalization such that a line search is used in the event that 
the step is initially rejected by the predicted versus actual decrease comparison. If the line 
search fails to find a viable step length for the Newton step, it falls back onto a scaled gradient 
or a gradient descent step. The trust radius is then modified based on the line search step length.

\subsection{Bounded Nonlinear Conjugate Gradient}
\label{sec:bncg}

BNCG extends the unconstrained nonlinear conjugate gradient algorithm to 
bound constraints via gradient projections and a bounded More-Thuente line search. 

Like its unconstrained counterpart, BNCG offers gradient descent and a variety of CG updates: 
Fletcher-Reeves, Polak-Ribi\'ere, Polak-Ribi\'ere-Plus, Hestenes-Stiefel, Dai-Yuan,
Hager-Zhang, Dai-Kou, Kou-Dai, and the Self-Scaling Memoryless (SSML) BFGS, DFP,
and Broyden methods.
These methods can be specified by using the command line argument 
{\tt -tao\_bncg\_type <gd,fr,pr,prp,hs,dy,hz,dk,kd,ssml\_bfgs,ssml\_dfp,ssml\_brdn>},
respectively. The default value is {\tt ssml\_bfgs}. We have scalar preconditioning for
these methods, and it is controlled by the flag {\tt tao\_bncg\_alpha}. To disable
rescaling, use $\alpha = -1.0$, otherwise $\alpha \in [0, 1]$.
BNCG is available via the TAO solver {\tt TAOBNCG} or the {\tt -tao\_type bncg} flag.

Some individual methods also contain their own parameters. The Hager-Zhang and Dou-Kai
 methods have a parameter that determines the minimum amount of contribution the previous
 search direction gives to the next search direction. The flags are
{\tt -tao\_bncg\_hz\_eta} and {\tt -tao\_bncg\_dk\_eta}, and by default are set to $0.4$ 
and $0.5$ respectively. The Kou-Dai method has multiple parameters.
{\tt -tao\_bncg\_zeta} serves the same purpose as the previous two; set to $0.1$ 
by default. There is also a parameter to scale the contribution of 
$y_k \equiv \nabla f(x_k) - \nabla f(x_{k-1}) $ in the search direction update. 
It is controlled by {\tt -tao\_bncg\_xi}, and is equal to $1.0$ by default. 
There are also times where we want to maximize the descent as measured 
by $\nabla f(x_k)^T d_k$, and that may be done by using a negative value of 
$\xi$; this achieves better performance when not using the diagonal preconditioner
described next. This is enabled by default, and is controlled by 
{\tt -tao\_bncg\_neg\_xi}. Finally, the Broyden method has its convex combination
parameter, set with {\tt -tao\_bncg\_theta}. We have this as 1.0 by 
default, i.e. it is by default the BFGS method. One can also individually 
tweak the BFGS and DFP contributions using the multiplicative constants 
{\tt -tao\_bncg\_<bfgs,dfp>\_scale}; both are set to $1$ by default.

All methods can be scaled using the parameter {\tt -tao\_bncg\_alpha}, which continuously
 varies in $[0, 1]$. The default value is set depending on the method from initial testing.

BNCG also offers a special type of method scaling. It employs Broyden diagonal scaling
as an option for its CG methods, turned on with the flag {\tt -tao\_bncg\_diag\_scaling}. Formulations for both the forward (regular) and inverse Broyden methods are developed,
 controlled by the flag {\tt -tao\_bncg\_mat\_lmvm\_forward}. It is set to True 
by default. Whether one uses the forward
or inverse formulations depends on the method being used. For example, in our preliminary computations,
the forward formulation works better for the SSML\_BFGS method, but the inverse formulation
works better for the Hestenes-Stiefel method. The convex combination parameter for the
 Broyden scaling is controlled by {\tt -tao\_bncg\_mat\_lmvm\_theta}, and is 0 by default. We also
employ rescaling of the Broyden diagonal, which aids the linesearch immensely. 
The rescaling parameter is controlled by {\tt -tao\_bncg\_mat\_lmvm\_alpha}, 
and should be $\in [0, 1]$. One can disable rescaling of the Broyden diagonal entirely by
 setting {\tt -tao\_bncg\_mat\_lmvm\_sigma\_hist 0}.

One can also supply their own preconditioner, serving as a Hessian initialization to
 the above diagonal scaling.
The appropriate user function in the code is {\tt TaoBNCGSetH0(tao, H0)} where 
{\tt H0} is the user-defined {\tt Mat} object that serves as a preconditioner. 
For an example of similar usage, see {\tt tao/examples/tutorials/ex3.c}.

The active set estimation uses the Bertsekas-based method described in Section~\ref{sec:bnk},
which can be deactivated using {\tt -tao\_bncg\_as\_type none}, in which case the 
algorithm will use the current iterate to determine the bounded variables with no tolerancing 
and no look-ahead step. As in the BNK algorithm, the initial bound tolerance and estimator
step length used in the Bertsekas method can be set via {\tt -tao\_bncg\_as\_tol} and 
{\tt -tao\_bncg\_as\_step}, respectively.

In addition to automatic scaled gradient descent restarts under certain local curvature conditions,
we also employ restarts based on a check on descent direction 
such that $\nabla f(x_k)^T d_k \in [-10^{11}, -1^{-9}]$.
Furthermore, we allow for a variety of alternative restart strategies, 
all disabled by default.
The {\tt -tao\_bncg\_unscaled\_restart} flag allows one to 
disable rescaling of the gradient for
gradient descent steps. The {\tt -tao\_bncg\_spaced\_restart} flag tells the 
solver to restart every $Mn$ iterations,
where $n$ is the problem dimension and $M$ is a constant 
determined by {\tt -tao\_bncg\_min\_restart\_num} and is
6 by default. We also have dynamic restart strategies based on 
checking if a function is locally quadratic; if so,
go do a gradient descent step. The flag is {\tt -tao\_bncg\_dynamic\_restart}, 
disabled by default since the CG
solver usually does better in those cases anyway. 
The minimum number of quadratic-like steps before a restart is 
set using {\tt -tao\_bncg\_min\_quad} and is 6 by default. 


\subsection{Trust-Region Newton Method}\label{sec:tron}\sindex{trust region}
The TRON \cite{lin_c3} algorithm is an active-set method that uses a 
combination of gradient projections and 
a preconditioned conjugate gradient method to minimize an objective function.
Each iteration of the TRON algorithm requires function, gradient, 
and Hessian evaluations.  In each iteration, the algorithm
first applies several conjugate gradient iterations.  
After these iterates, the TRON solver momentarily ignores the variables
that equal one of its bounds and
applies a preconditioned conjugate gradient method to a
quadratic model of the remaining set of ``free'' variables.  

The TRON algorithm solves a reduced linear system
defined by the rows and columns corresponding to the variables that
lie between the upper and lower bounds.
The TRON algorithm applies a trust region to the 
conjugate gradients to ensure convergence.  The initial trust-region
radius
can be set by using the command 
{\tt Tao\-Set\-Initial\-Trust\-Region\-Radius()},
and the current trust region size can be found by using the command
{\tt Tao\-Get\-Current\-Trust\-Region\-Radius()}.
The initial trust region can significantly alter the 
rate of convergence for the algorithm and should be
tuned and adjusted for optimal performance.

\subsection{Bound-constrained Limited-Memory Variable-Metric Method}
\label{sec:blmvm}

BLMVM is a limited-memory, variable-metric method and is 
the bound-constrained variant of the LMVM method for
unconstrained optimization. It uses projected gradients to approximate
the Hessian, eliminating the need for Hessian evaluations.
The method can be set by using the TAO solver {\tt tao\_blmvm}. For more details, 
please see the LMVM section in the unconstrained algorithms as well as the 
LMVM matrix documentation in the PETSc manual.

This algorithm will be deprecated in the next version in favor of the bounded quasi-Newton 
line search (BQNLS) algorithm.

\subsection{Bounded Quasi-Newton-Krylov}\label{sec:bqnk}

BQNK algorithms use the BNK infrastructure, but replace the exact Hessian with a quasi-Newton 
approximation. The matrix-free forward product operation based on quasi-Newton update formulas 
are used in conjunction with Krylov solvers to compute step directions. The quasi-Newton inverse 
application is used to precondition the Krylov solution, and typically helps converge to a step 
direction in $\mathcal{O}(10)$ iterations. This approach is most useful with quasi-Newton update 
types such as Symmetric Rank-1 that cannot strictly guarantee positive-definiteness. The BNLS 
framework with Hessian shifting, or the BNTR framework with trust region safeguards, can 
successfully compensate for the Hessian approximation becoming indefinite.

Similar to the full Newton-Krylov counterpart, BQNK algorithms come in three forms separated by 
the globalization technique: line search (BQNKLS), trust region (BQNKTR) and trust region w/ line 
search fall-back (BQNKTL). These algorithms are available via {\tt tao\_type <bqnkls, bqnktr, bqnktl>}.

\subsection{Bounded Quasi-Newton Line Search (BQNLS)}\label{sec:bqnls}

BQNLS algorithm uses the BNLS infrastructure, but replaces the step calculation with a direct 
inverse application of the approximate Hessian based on quasi-Newton update formulas. No Krylov 
solver is used in the solution, and therefore the quasi-Newton method chosen must guarantee 
a positive-definite Hessian approximation. This algorithm is available via {\tt tao\_type bqnls}.

\section{PDE-Constrained Optimization}
\label{sec:lcl}

TAO solves PDE-constrained optimization problems of the form
\[
\begin{array}{ll}
\displaystyle \min_{u,v} & f(u,v) \\
\mbox{subject to} & g(u,v) = 0,
\end{array}
\]
where the state variable $u$ is the solution to the discretized
partial differential equation defined by $g$ and parametrized by 
the design variable $v$, and $f$ is an objective function.  The 
Lagrange multipliers on the constraint are denoted by $y$.
This method is set by using the linearly constrained augmented
Lagrangian TAO solver {\tt tao\_lcl}.

We make two main assumptions when solving these problems: the 
objective function and PDE constraints 
have been discretized so that we can treat the optimization problem 
as finite dimensional and $\nabla_u g(u,v)$ is invertible for all 
$u$ and $v$.
% If required, we can assume $\|\nabla_u g(u,v)\|$ is uniformly bounded 
% above and below in an appropriate norm.

\subsection{Linearly-Constrained Augmented Lagrangian Method}
Given the current iterate $(u_k, v_k, y_k)$, the linearly constrained
augmented Lagrangian method approximately solves the optimization 
problem
\[
\begin{array}{ll}
\displaystyle \min_{u,v} & \tilde{f}_k(u, v) \\
\mbox{subject to} & A_k (u-u_k) + B_k (v-v_k) + g_k = 0,
\end{array}
\]
where $A_k = \nabla_u g(u_k,v_k)$, $B_k = \nabla_v g(u_k,v_k)$, and
$g_k = g(u_k, v_k)$ and 
\[
\tilde{f}_k(u,v) = f(u,v) - g(u,v)^T y^k + \frac{\rho_k}{2} \| g(u,v) \|^2
\]
is the augmented Lagrangian function.  This optimization problem is
solved in two stages.  The first computes the Newton direction and
finds a feasible point for the linear constraints.  The second 
computes a reduced-space direction that maintains feasibility
with respect to the linearized constraints and improves the
augmented Lagrangian merit function.

\subsubsection{Newton Step}

The Newton direction is obtained by fixing the design variables at
their current value and solving the linearized constraint for the
state variables.  In particular, we solve the system of equations
\[
A_k du = -g_k
\]
to obtain a direction $du$.  We need a direction that provides
sufficient descent for the merit function
\[
  \frac{1}{2} \|g(u,v)\|^2.
\]
That is, we require $g_k^T A_k du < 0$.  

If the Newton direction is a descent direction, then we choose a 
penalty parameter $\rho_k$ so that $du$ is also a sufficient descent 
direction for the augmented Lagrangian merit function.  We 
then find $\alpha$ to approximately minimize the augmented 
Lagrangian merit function along the Newton direction.
\[
\displaystyle \min_{\alpha \geq 0} \; \tilde{f}_k(u_k + \alpha du, v_k).
\]
We can enforce either the sufficient decrease condition or the 
Wolfe conditions during the search procedure.  The new point, 
\[
\begin{array}{lcl}
u_{k+\frac{1}{2}} & = & u_k + \alpha_k du \\
v_{k+\frac{1}{2}} & = & v_k,
\end{array}
\]
satisfies the linear constraint
\[
A_k (u_{k+\frac{1}{2}} - u_k) + B_k (v_{k+\frac{1}{2}} - v_k) + \alpha_k g_k = 0.
\]
If the Newton direction computed does not provide descent for the merit 
function, then we can use the steepest descent direction $du = -A_k^T g_k$ 
during the search procedure.  However, the implication that the intermediate
point approximately satisfies the linear constraint is no longer true.
% We many need to apply a feasibility restoration phase until the Newton
% direction provided sufficient descent for the merit function.
% Note that we can update a limited memory quasi-Newton approximation
% to $\nabla^2_{u,u} \tilde{f}_k$ as this point, but doing so does not
% seem to be useful, unless it is being used perhaps to correct the 
% multiplier estimates.

\subsubsection{Modified Reduced-Space Step}
We are now ready to compute a reduced-space step for the modified
optimization problem:
\[
\begin{array}{ll}
\displaystyle \min_{u,v} & \tilde{f}_k(u, v) \\
\mbox{subject to} & A_k (u-u_k) + B_k (v-v_k) + \alpha_k g_k = 0.
\end{array}
\]
We begin with the change of variables
\[
\begin{array}{ll}
\displaystyle \min_{du,dv} & \tilde{f}_k(u_k+du, v_k+dv) \\
\mbox{subject to} & A_k du + B_k dv + \alpha_k g_k = 0
\end{array}
\]
and make the substitution
\[
  du = -A_k^{-1}(B_k dv + \alpha_k g_k).
\]
Hence, the unconstrained optimization problem we need to solve is
\[
\begin{array}{ll}
\displaystyle \min_{dv} & \tilde{f}_k(u_k-A_k^{-1}(B_k dv + \alpha_k g_k), v_k+dv), \\
\end{array}
\]
which is equivalent to
\[
\begin{array}{ll}
\displaystyle \min_{dv} & \tilde{f}_k(u_{k+\frac{1}{2}} - A_k^{-1} B_k dv, v_{k+\frac{1}{2}}+dv). \\
\end{array}
\]
We apply one step of a limited-memory quasi-Newton method to this
problem.  The direction is obtain by solving the quadratic problem
\[
\begin{array}{ll}
\displaystyle \min_{dv} & \frac{1}{2} dv^T \tilde{H}_k dv + \tilde{g}_{k+\frac{1}{2}}^T dv,
\end{array}
\]
where $\tilde{H}_k$ is the limited-memory quasi-Newton approximation
to the reduced Hessian matrix, a positive-definite matrix, and 
$\tilde{g}_{k+\frac{1}{2}}$ is the reduced gradient.
\[
\begin{array}{lcl}
\tilde{g}_{k+\frac{1}{2}} & = & \nabla_v \tilde{f}_k(u_{k+\frac{1}{2}}, v_{k+\frac{1}{2}}) -
	      \nabla_u \tilde{f}_k(u_{k+\frac{1}{2}}, v_{k+\frac{1}{2}}) A_k^{-1} B_k \\
	   & = & d_{k+\frac{1}{2}} + c_{k+\frac{1}{2}} A_k^{-1} B_k
\end{array}
\]
The reduced gradient is obtained from one linearized adjoint solve
\[
y_{k+\frac{1}{2}} = A_k^{-T}c_{k+\frac{1}{2}}
\]
and some linear algebra
\[
\tilde{g}_{k+\frac{1}{2}} = d_{k+\frac{1}{2}} + y_{k+\frac{1}{2}}^T B_k.
\]
Because the Hessian approximation is positive definite and we know its
inverse, we obtain the direction
\[
  dv = -H_k^{-1} \tilde{g}_{k+\frac{1}{2}}
\]
and recover the full-space direction from one linearized forward solve,
\[
  du = -A_k^{-1} B_k dv.
\]
Having the full-space direction, which satisfies the linear constraint, 
we now approximately minimize the augmented Lagrangian merit function 
along the direction.
\[
\begin{array}{lcl}
\displaystyle \min_{\beta \geq 0} & \tilde{f_k}(u_{k+\frac{1}{2}} + \beta du, v_{k+\frac{1}{2}} + \beta dv)
\end{array}
\]
We enforce the Wolfe conditions during the search procedure.  The new point
is
\[
\begin{array}{lcl}
u_{k+1} & = & u_{k+\frac{1}{2}} + \beta_k du \\
v_{k+1} & = & v_{k+\frac{1}{2}} + \beta_k dv.
\end{array}
\]
The reduced gradient at the new point is computed from
\[
\begin{array}{lcl}
y_{k+1} & = & A_k^{-T}c_{k+1} \\
\tilde{g}_{k+1} & = & d_{k+1} - y_{k+1}^T B_k,
\end{array}
\]
where $c_{k+1} = \nabla_u \tilde{f}_k (u_{k+1},v_{k+1})$ and
$d_{k+1} = \nabla_v \tilde{f}_k (u_{k+1},v_{k+1})$.  The
multipliers $y_{k+1}$ become the multipliers used in the
next iteration of the code.  The quantities $v_{k+\frac{1}{2}}$,
$v_{k+1}$, $\tilde{g}_{k+\frac{1}{2}}$, and $\tilde{g}_{k+1}$ are
used to update $H_k$ to obtain the limited-memory quasi-Newton
approximation to the reduced Hessian matrix used in the next
iteration of the code.  The update is skipped if it cannot be
performed.

\section{Nonlinear Least-Squares}
\label{sec:leastsquares}
Given a function $F: \Re^n \to \Re^m$, the nonlinear least-squares 
problem minimizes 
\begin{equation} 
 f(x)= \| F(x) \|_2^2 = \sum_{i=1}^m F_i(x)^2.
 \label{eq:nlsf}
\end{equation}
The nonlinear equations $F$ should be specified with the function
{\tt Tao\-Set\-Separ\-able\-Objective\-Function()}.
Although the current version of TAO has no derivative-based algorithms 
for nonlinear least squares, we anticipate that such an algorithm will 
be available in the next release. 

\begin{comment}
Most solvers also require the Jacobian of this function.  
The routine that evaluates this matrix should be set
using the command {\tt Tao\-Set\-Jacobian\-Routine()}.

Using the constraint data provided by the user, the nonlinear
least squares solver formulates a minimization problem, and
solves the minimization problem using another TAO solver.  
A nonlinear least squares application should be solved with
the method \texttt{tao\_pounders} or \texttt{tao\_lm}.
\textbf{***** Jason, check the lm above. It would be great to include it
but otherwise say that we only have methods that do not use the Jacobian
and that if they have a Jacobian, they should use something else/the
next/beta version.}
\end{comment}

\subsection{POUNDerS}\label{sec:pounders}
One algorithm for solving the least squares problem (\ref{eq:nlsf}) when
the Jacobian of the residual vector $F$ is unavailable is the model-based
POUNDerS (Practical Optimization Using No Derivatives for
sums of Squares) algorithm (\texttt{tao\_pounders}). 
POUNDerS employs a derivative-free trust-region framework as described in
\cite{Dfobook} in order to converge to local minimizers. An example of
this version of POUNDerS applied to a
practical least-squares problem can be found in \cite{UNEDF0}.

\subsubsection{Derivative-Free Trust-Region Algorithm}
In each iteration $k$, the algorithm maintains a model $m_k(x)$,
described below, of the nonlinear least squares function $f$ centered about
the current iterate $x_k$. 

If one assumes that the maximum number of function evaluations has not been
reached and that $\|\nabla m_k(x_k)\|_2>$\texttt{gtol}, the next point 
$x_+$ to be evaluated is obtained by solving the trust-region subproblem 
\begin{equation}
 \min\left\{
 m_k(x) : 
 \|x-x_k\|_{p} \leq \Delta_k, 
 \right \}, 
 \label{eq:poundersp}
\end{equation}
where $\Delta_k$ is the current 
trust-region radius. By default we use a trust-region norm with $p=\infty$
and solve (\ref{eq:poundersp}) with the BLMVM method described in
Section~\ref{sec:blmvm}.  While the subproblem is a bound-constrained
quadratic program, it may not be convex and the BQPIP and GPCG methods 
may not solve the subproblem.  Therefore, either BLMVM (the default) 
or TRON should be used.  Note: TRON uses its own internal trust region 
that may interfere with the infinity-norm trust region used in the 
model problem (\ref{eq:poundersp}).

The residual vector is then evaluated to obtain $F(x_+)$ and hence
$f(x_+)$. The ratio of actual decrease to predicted decrease, 
\[
 \rho_k = \frac{f(x_k)-f(x_+)}{m_k(x_k)-m_k(x_+)},
\]
as well as an indicator, \texttt{valid}, on the model's quality of
approximation on the trust region is then used to update the iterate,
\[
x_{k+1} = \left\{\begin{array}{ll}
x_+ & \mbox{if } \rho_k \geq \eta_1 \\
x_+ & \mbox{if } 0<\rho_k <\eta_1  \mbox{ and \texttt{valid}=\texttt{true}}
\\
x_k & \mbox{else}, 
\end{array}
\right.
\]
and trust-region radius,
\[
\Delta_{k+1} = \left\{\begin{array}{ll}
 \mbox{min}(\gamma_1\Delta_k, \Delta_{\max}) & \mbox{if } \rho_k \geq
\eta_1 \mbox{ and } \|x_+-x_k\|_p\geq \omega_1\Delta_k \\
\gamma_0\Delta_k & \mbox{if } \rho_k < \eta_1 \mbox{ and
\texttt{valid}=\texttt{true}} \\
\Delta_k &  \mbox{else,}
\end{array}
\right.
\]
where $0 < \eta_1 < 1$,  $0 < \gamma_0 < 1 < \gamma_1$, $0<\omega_1<1$,
and $\Delta_{\max}$ are constants. 

If $\rho_k\leq 0$ and \texttt{valid} is \texttt{false}, the
iterate and trust-region radius remain unchanged after the above updates,
and the algorithm tests whether the direction $x_+-x_k$ improves the model. If
not, the algorithm performs an additional evaluation to obtain
$F(x_k+d_k)$, where $d_k$ is a model-improving direction. 

The iteration counter is then updated, and the next model $m_{k}$ is
obtained as described next.

\subsubsection{Forming the Trust-Region Model}
In each iteration, POUNDerS uses a subset of the available evaluated
residual vectors $\{ F(y_1), F(y_2), \cdots \}$ to form an interpolatory
quadratic model of each residual component. The $m$ quadratic models
\begin{equation}
 q_k^{(i)}(x) = 
 F_i(x_k) + (x-x_k)^T g_k^{(i)} + \frac{1}{2} (x-x_k)^T H_k^{(i)} (x-x_k), 
 \qquad i = 1, \ldots, m
\label{eq:models}
\end{equation}
thus satisfy the interpolation conditions 
\[ 
 q_k^{(i)}(y_j) = F_i(y_j), \qquad i=1, \ldots, m; \, j=1,\ldots , l_k
\]
on a common interpolation set $\{y_1, \cdots , y_{l_k}\}$ of size
$l_k\in[n+1,$\texttt{npmax}$]$.

The gradients and Hessians of the models in (\ref{eq:models}) are then
used to construct the master model, 
\begin{equation}
m_k(x) = f(x_k) + 
2(x-x_k)^T \sum_{i=1}^{m} F_i(x_k) g_k^{(i)} +
(x-x_k)^T \sum_{i=1}^{m}
\left( g_k^{(i)} \left(g_k^{(i)}\right)^T +  F_i(x_k)
H_k^{(i)}\right) (x-x_k).
 \label{eq:newton2}
\end{equation}
The process of forming these models also computes the indicator
\texttt{valid} of the model's local quality.


\subsubsection{Parameters}
\begin{comment}
If evaluations of the residual vector $F$ are available, then this can be
input and will often result in improved performance. Assuming that
\texttt{nfs}$>0$ such points are available, then the user should input the
matrix of locations \texttt{X0} $=[y_1, \ldots, y_{\rm nfs}]$ and the
corresponding matrix of residuals \texttt{F0} $= [F(y_1), \ldots, F(y_{\rm
nfs})]$. The user selects which one of these points is the starting point
(this is typically the point $\arg \min_{i=1,\ldots, {\rm nfs}}
\|F(y_i)\|^2$) to be used and inputs the index as the integer
\texttt{xkin}.

If no such points are available, then \texttt{nfs}$=0$, \texttt{F0} is
empty, \texttt{xkin}$=1$, and the user inputs the desired starting point
(whose residual has not yet been evaluated) as \texttt{X0} $=[x_0]$.


In addition to the inputs \texttt{nfs}, \texttt{X0}, \texttt{F0},
\texttt{xkin} above, 
\end{comment}

POUNDerS supports the following parameters that can be set from the
command line or PETSc options file:
\begin{description}
  \item[\texttt{-tao\_pounders\_delta <delta>}] The initial trust-region 
  radius ($>0$, real).
  This is  used to determine the size of the initial neighborhood within
  which the algorithm should look. 
  %In particular, if \texttt{nfs}=0, the
  %algorithm will evaluate $F$ at $n$ points within a distance
  %\texttt{delta} of \texttt{X0};
  
  \item[\texttt{-tao\_pounders\_npmax <npmax>}] The maximum number of 
  interpolation points used
  ($n+2\leq$ \texttt{npmax} $\leq 0.5(n+1)(n+2)$). This input is made
  available to advanced users. We recommend the default
  value (\texttt{npmax}$=2n+1$) be used by others.

  \item[\texttt{-tao\_pounders\_gqt}] Use the gqt algorithm to solve the
  subproblem~(\ref{eq:poundersp}) (uses $p=2$) instead of BQPIP.

  \item[\texttt{-pounders\_subsolver}] If the default
BQPIP algorithm is used to solve the 
  subproblem~(\ref{eq:poundersp}), the parameters of the subproblem solver
  can be accessed using the command line options  prefix
\texttt{-pounders\_subsolver\_}.  For example, 
\begin{verbatim}
       -pounders_subsolver_tao_gatol 1.0e-5
\end{verbatim}  
sets the gradient tolerance of the subproblem solver to $10^{-5}$.
\end{description}

Additionally,  the user provides an initial solution vector, a vector for
storing
the separable objective function, and a routine for evaluating the residual
vector $F$.  These are described in detail in 
Sections~\ref{sec:fghj} and \ref{sec:evalsof}. 
Here we remark that because gradient information is not available for
scaling purposes, it can be useful to ensure that the problem is
reasonably well scaled. A simple way to do so is to rescale the
decision variables $x$ so that their typical values are expected to lie
within the unit hypercube $[0,1]^n$.

\subsubsection{Convergence Notes}
Because the gradient function is not provided to POUNDerS, 
the norm of the gradient of the objective function is not available.
Therefore, 
for convergence criteria, this norm is approximated by the norm of the
model gradient and used only when the model gradient is deemed to be a
reasonable approximation of the gradient of the objective. In practice, the
typical grounds for termination for expensive derivative-free problems is
the maximum number of function evaluations allowed. 




\section{Complementarity}
\label{sec:complementarity}

Mixed complementarity problems, or box-constrained variational inequalities,
are related to nonlinear systems of equations.  They are defined by a 
continuously differentiable function,
$F:\Re^n \to \Re^n$, and bounds, $\ell \in \{\Re \cup \{-\infty\}\}^n$ 
and $u \in \{\Re \cup \{\infty\}\}^n$, on the variables such that 
$\ell \leq u$.  Given this information, $\x^* \in [\ell,u]$ is a solution to 
MCP($F$, $\ell$, $u$) if for each $i \in \{1, \ldots, n\}$ we have at
least one of the following:
\begin{eqnarray*}
\begin{array}{ll}
F_i(x^*) \geq 0 & \mbox{if } x^*_i = \ell_i \\
F_i(x^*) = 0 & \mbox{if } \ell_i < x^*_i < u_i \\
F_i(x^*) \leq 0 & \mbox{if } x^*_i = u_i.
\end{array}
\end{eqnarray*}
Note that when $\ell = \{-\infty\}^n$ and $u = \{\infty\}^n$, we have a 
nonlinear system of equations, and $\ell = \{0\}^n$ and $u = \{\infty\}^n$ 
correspond to the nonlinear complementarity problem \cite{cottle:nonlinear}.

Simple complementarity conditions arise from the first-order optimality 
conditions from optimization \cite{karush:minima, kuhn.tucker:nonlinear}.  
In the simple bound-constrained optimization case, these conditions 
correspond to MCP($\nabla f$, $\ell$, $u$), where $f: \Re^n \to \Re$ 
is the objective function.  In a one-dimensional setting these conditions 
are intuitive.  If the solution is at the lower bound, then the function must 
be increasing and $\nabla f \geq 0$.  If the solution is at the 
upper bound, then the function must be decreasing and $\nabla f \leq 0$.  
If the solution 
is strictly between the bounds, we must be at a stationary point and 
$\nabla f = 0$.  Other complementarity problems arise in economics and 
engineering \cite{ferris.pang:engineering}, game 
theory \cite{nash:equilibrium}, and finance \cite{huang.pang:option}.

Evaluation routines for $F$ and its Jacobian must be supplied prior
to solving the application.
The bounds, $[\ell,u]$, on the variables must also be 
provided.  
If no starting point is supplied, a default starting point of all zeros 
is used.

\subsection{Semismooth Methods}

TAO has two implementations of semismooth algorithms 
\cite{munson.facchinei.ea:semismooth, deluca.facchinei.ea:semismooth, 
facchinei.fischer.ea:semismooth} for solving mixed complementarity problems.
Both are based on a reformulation of the mixed complementarity problem
as a nonsmooth system of equations using the Fischer-Burmeister 
function \cite{fischer:special}.  A nonsmooth Newton method is applied to
the reformulated system to calculate a solution.  The theoretical properties
of such methods are detailed in the aforementioned references.

The Fischer-Burmeister function, $\phi:\Re^2 \to \Re$, is defined as
\begin{eqnarray*}
\phi(a,b) := \sqrt{a^2 + b^2} - a - b.
\end{eqnarray*}
This function has the following key property,
\begin{eqnarray*}
\begin{array}{lcr}
        \phi(a,b) = 0 & \Leftrightarrow & a \geq 0,\; b \geq 0,\; ab = 0,
\end{array}
\end{eqnarray*}
used when reformulating the mixed complementarity problem as the system of 
equations $\Phi(x) = 0$, where $\Phi:\Re^n \to \Re^n$.  
The reformulation is defined componentwise as
\begin{eqnarray*}
\Phi_i(x) := \left\{ \begin{array}{ll}
   \phi(x_i - l_i, F_i(x)) & \mbox{if } -\infty < l_i < u_i = \infty, \\
   -\phi(u_i-x_i, -F_i(x)) & \mbox{if } -\infty = l_i < u_i < \infty, \\
   \phi(x_i - l_i, \phi(u_i - x_i, - F_i(x))) & \mbox{if } -\infty < l_i < u_i < \infty, \\
   -F_i(x) & \mbox{if } -\infty = l_i < u_i = \infty, \\
   l_i - x_i & \mbox{if } -\infty < l_i = u_i < \infty.
   \end{array} \right.
\end{eqnarray*}
We note that $\Phi$ is not differentiable everywhere but satisfies a
semismoothness property 
\cite{mifflin:semismooth, qi:convergence, qi.sun:nonsmooth}.  Furthermore,
the natural merit function, $\Psi(x) := \frac{1}{2} \| \Phi(x) \|_2^2$, is 
continuously differentiable.

The two semismooth TAO solvers both solve the system $\Phi(x) = 0$ by applying
a nonsmooth Newton method with a line search.  We calculate a direction, $d^k$,
by solving the system $H^kd^k = -\Phi(x^k)$, where $H^k$ is an element of the 
$B$-subdifferential \cite{qi.sun:nonsmooth} of $\Phi$ at $x^k$.  If the 
direction calculated does not satisfy a suitable descent condition, then 
we use the negative gradient of the merit function, $-\nabla \Psi(x^k)$, as 
the search direction.  A standard Armijo search \cite{armijo:minimization} is
used to find the new iteration.  Nonmonotone searches 
\cite{grippo.lampariello.ea:nonmonotone} are also available by setting 
appropriate runtime options.  See Section \ref{sec:TaoLineSearch} for further 
details.

The first semismooth algorithm available in TAO is not guaranteed to 
remain feasible with respect to the bounds, $[\ell, u]$, and is termed 
an infeasible semismooth method.  This method can be specified by using the 
{\tt tao\_ssils} solver.  In this case, the descent test used is 
that
\begin{eqnarray*}
\nabla \Psi(x^k)^Td^k \leq -\delta\| d^k \|^\rho.
\end{eqnarray*}
Both $\delta > 0$ and $\rho > 2$ can be modified by using the
runtime options {\tt -tao\_ssils\_delta <delta>} and 
{\tt -tao\_ssils\_rho <rho>}, respectively. By default, 
$\delta = 10^{-10}$ and $\rho = 2.1$.

An alternative is to remain feasible with respect to the bounds by using a 
projected Armijo line search.  This method can be specified by using the 
{\tt tao\_ssfls} solver.  The descent test used is the same as above 
where the direction in this case corresponds to the first part of the 
piecewise linear arc searched by the projected line search.  
Both $\delta > 0$ and $\rho > 2$ can be modified by using the runtime 
options  {\tt -tao\_ssfls\_delta <delta>} and {\tt -tao\_ssfls\_rho <rho>} 
respectively.  By default, $\delta = 10^{-10}$ and $\rho = 2.1$.

The recommended algorithm is the infeasible semismooth method,
{\tt tao\_ssils}, because of its strong global and local convergence 
properties.  However, if it is known that $F$ is not defined outside 
of the box, $[\ell,u]$, perhaps because of the presence of $\log$ functions,
the feasibility-enforcing version of the algorithm, {\tt tao\_ssfls}, is a 
reasonable alternative.

\section{Quadratic Solvers}

Quadratic solvers solve optimization problems of the form
\[
\begin{array}{ll}
\displaystyle \min_{x} & \frac{1}{2}x^T Q x + c^T x \\
\mbox{subject to} & l \geq x \geq u
\end{array}
\]
where the gradient and the Hessian of the objective are both constant.

\subsection{Gradient Projection Conjugate Gradient Method}
The GPCG \cite{more-toraldo} algorithm is much like the TRON algorithm, discussed in
Section \ref{sec:tron}, except that
it assumes that the objective function is quadratic and convex.
Therefore, it evaluates the function, gradient, and Hessian only
once.
Since the objective function
is quadratic, the algorithm does not use a trust region.  
All the options that apply to TRON except for trust-region
options also apply to GPCG. It can be set by using the TAO solver 
{\tt tao\_gpcg} or via the optio flag {\tt -tao\_type gpcg}.

\subsection{Interior-Point Newton's Method}\label{sec:bqpip}
The BQPIP algorithm is an interior-point method for bound
constrained quadratic optimization.  It can be set by using the
TAO solver of {\tt tao\_bqpip} or via the option flag {\tt -tao\_type bgpip}.
Since it assumes the objective function is quadratic, 
it evaluates the function, gradient, and Hessian only once.
This method also requires the solution of systems of linear equations,
whose solver can be accessed and modified 
with the command {\tt Tao\-Get\-KSP()}.
